#include "includes2DPC.h"
#include "pocs2DPC.h"
// Code performing Hilbert projection.
using namespace std;
void P1hilbert(ScanData & s, int j)
{
	int i,n;
	float pterm=0,x=0, temp =0, temp2=0,a = 0.5 ,b = 1.-a;
	vector<float> image (s.Ny, 0);
	
	// Compute preterm
	for (i=0; i<s.Nx;i++){
		if(s.CSmask[i][j]) pterm = pterm + s.dx*s.recon[i][j];
	}
	//cout<< pterm<<endl;
	// Compute H transform


	for(i=0; i<s.Nx;++i){ // x loop
		image[i]=0;
		for(n=0; n<s.Nx;++n){ // x prime loop
			if( (abs(i-n)%2 != 0) ) image[i] = image[i] + s.recon[n][j] * (2./(n-i))/M_PI;
		}
		s.HT[i][j] = image[i];
	}
	// perform smooting
	temp = image[0];
	for(i=1; i<s.Nx-1;++i){ // x loop
		temp2=image[i];
		image[i] = 0.5*b*temp + a*image[i] + 0.5*b*image[i+1];
		temp=temp2;
	}

	// Compute projection ||((i-n) !=0)
	for(i=0; i<s.Nx;++i){
		if (s.DBPmask[i][j]) image[i] = s.DBPdata[i][j]*(-0.5/M_PI) + 0.019;
	}
	
	// Compute inverse H transform
	for(i=0; i<s.Nx;++i){ // preweight
		x = s.xmin + i*s.dx;
		image[i] = image[i] * sqrt(s.X2 - x*x);
	}
	
	for(i=0; i<s.Nx;++i){ // convolution
		s.recon[i][j]=0;
		for(n=0; n<s.Nx;++n){ // x prime loop
			if( abs(i-n)%2 != 0) s.recon[i][j] = s.recon[i][j] + image[n] * (-2./(n-i));
		}
		
	}

	// perform smooting
	temp = image[0];
	for(i=1; i<s.Nx-1;++i){ // x loop
		temp2=image[i];
		image[i] = 0.25*temp + 0.5*image[i] + 0.25*image[i+1];
		temp=temp2;
	}

	for(i=0; i<s.Nx;++i){ // postweight
		x = s.xmin + i*s.dx;
		s.recon[i][j] = (pterm + s.recon[i][j])/sqrt(s.X2 - x*x)/M_PI;
		//s.recon[i][j] = (s.recon[i][j])/sqrt(s.X2 - x*x)/M_PI;
		image[i]=0;
	}
	
	pterm =0;
}


void P1hilbert_perp(ScanData & s, int i)
{
	int j,n;
	float pterm=0,y=0, temp =0, temp2=0,a = 0.50 ,b = 1.-a;
	vector<float> image (s.Ny, 0);
	
	// Compute preterm
	for (j=0; j<s.Ny;j++){
		if(s.CSmask[i][j]) pterm = pterm + s.dy*s.recon[i][j];
	}
	//cout<< pterm<<endl;
	// Compute H transform


	for(j=0; j<s.Ny;++j){ // y loop
		image[j]=0;
		for(n=0; n<s.Ny;++n){ // x prime loop
			if( (abs(j-n)%2 != 0) ) image[j] = image[j] + s.recon[i][n] * (2./(n-j))/M_PI;
		}
		s.HT[i][j] = image[j];
	}
	// perform smooting
	temp = image[0];
	for(j=1; j<s.Ny-1;++j){ // x loop
		temp2=image[j];
		image[j] = 0.5*b*temp + a*image[j] + 0.5*b*image[j+1];
		temp=temp2;
	}

	// Compute projection ||((i-n) !=0)
	for(j=0; j<s.Ny;++j){
		if (s.DBPmask[i][j]) image[j] = s.DBPdataPerp[i][j]*(0.5/M_PI) - 0.032;
	}
	
	// Compute inverse H transform
	for(j=0; j<s.Ny;++j){ // preweight
		y = s.ymin + j*s.dy;
		image[j] = image[j] * sqrt(s.X2 - y*y);
	}
	
	for(j=0; j<s.Ny;++j){ // convolution
		s.recon[i][j]=0;
		for(n=0; n<s.Ny;++n){ // x prime loop
			if( abs(j-n)%2 != 0) s.recon[i][j] = s.recon[i][j] + image[n] * (-2./(n-j));
		}
		
	}

	// perform smooting
	temp = image[0];
	for(j=1; j<s.Ny-1;++j){ // x loop
		temp2=image[j];
		image[j] = 0.25*temp + 0.5*image[j] + 0.25*image[j+1];
		temp=temp2;
	}

	for(j=0; j<s.Ny;++j){ // postweight
		y = s.ymin + j*s.dy;
		s.recon[i][j] = (pterm + s.recon[i][j])/sqrt(s.X2 - y*y)/M_PI;
		//s.recon[i][j] = (s.recon[i][j])/sqrt(s.X2 - x*x)/M_PI;
		image[j]=0;
	}
	
	pterm =0;
}
